*-----------------------------------------

// inputs are 
//
// alpha, lnTFP variable  and level of aggregation(macro input)
// l x (log of labor and land)

 cap program drop CalculateEfficiency
 program def CalculateEfficiency
  
  
  

	cap preserve
	 display "Calculate effiency gains using  $TFP as TFP measure,  aggregated at $aggregation level"
	
    cap drop parish_code
	quietly egen parish_code=group(region parish_name)
    quietly  drop if  parish_code==.
 
 	quietly  global gamma = $alphal+  $alphax 
 
	*5.2. Obtain levels of land and labor and tfp from logs (note all at plot level!) and totals (I calculate totals by year-season-village (parishID) )
	cap drop land
	cap drop labor
	quietly  gen land=exp(l)
	quietly  gen labor=exp(x)
	
 	quietly  cap drop TFP
	quietly  gen   TFP =exp( $lnTFP )
 
	
	quietly  bys   $aggregation: egen total_land=sum(land)
	quietly  bys $aggregation: egen total_labor=sum(labor)
 
	quietly gen z = TFP ^(1/(1-$gamma )) /100000000000000000000000  /*I do this transf. to reduce size of variable! does not affect shares... */
 
	 quietly bys $aggregation: egen total_z  =sum(z)
 
	*  obtain efficeicn allcoations and efficeint output  
 
	cap drop  land_effic labor_effic output_pred output_effic Y_pred Y_max
	
			quietly    gen land_effic =  total_land * z  / total_z 
 
			quietly    gen labor_effic =  total_labor * z  / total_z 
 			
			quietly  gen output_pred  =TFP  * ( land^$alphal  )*(labor^$alphax  )  
			
			quietly   gen output_effic =  TFP * (land_effic^$alphal  ) * (labor_effic ^$alphax  )
		
// 		    quietly  bys $aggregation: egen  Y_pred =sum(output_pred )
// 		   quietly bys  $aggregation: egen  Y_max =sum(output_effic )

  		    quietly  bys $time: egen  Y_pred =sum(output_pred )
 		   quietly bys  $time: egen  Y_max =sum(output_effic )

	quietly  gen effGain = (Y_max / Y_pred -1)*100
	bys $time: gen tmp_n=_n
	  


	  sum effGain if tmp_n==1
	  
	  
 	 quietly  sum $lnTFP, de
 display "Variance of ln TFP is " `r(Var)'
// 	 quietly  sum y, de
// 	 display "Stdev of ln output is " `r(sd)'
// 	 quietly  sum l, de
// 	 display "Stdev of ln land is " `r(sd)'
// 	 quietly  sum x, de
// 	 display "Stdev of ln labor is " `r(sd)'
//
	 	 quietly reghdfe l $lnTFP, abs( $time)
	 di "land -p roductivity elasticity is " _b[$lnTFP ]
		 quietly reghdfe  x $lnTFP, abs( $time)
	 di "labor - productivity elasticity is " _b[$lnTFP ]
	 	  quietly	  sum effGain
	di "No. production units ="	   `r(N)'
//	  
  cap restore
end

******************************************************************************************
******************************************************************************************



******************************************************************************************
******************************************************************************************


 cap program drop CalculateEfficiencyHH
 program def CalculateEfficiencyHH
  
  
  

	preserve
	
	*first aggregate data at HH levels
	
	gen land=exp(l)
	gen labor=exp(x)
		gen output=exp(y)

	  bys HHID year season: egen L_hh=sum(land)
	bys HHID year season: egen X_hh=sum(labor)
		bys HHID year season: egen Y_hh=sum(output)

//  	 *take weighted average of TFPA
// 	  gen phi_l=  land /  L_hh
// 	  gen phi_x=   labor / X_hh
// 
//   	quietly  gen   TFP =exp( $lnTFP )
//
//
// 		cap drop tmp
// 		gen tmp = TFP  *   phi_l^$alphal  * phi_x^$alphax 
// 	   bys HHID year season: egen s = sum(tmp)
// 	   cap drop $lnTFP 
// 	   	gen $lnTFP  =ln(s )

	cap drop $lnTFP 
 	gen $lnTFP  =ln(Y_hh ) -  $alphal * ln(L_hh) -  $alphax * ln(X_hh)
		
  *create identifier of HH sample
     bys HHID year season: gen n=_n
	 keep if  n==1
		drop x l y
	 gen l=ln(L_hh)
	 gen x=ln(X_hh)
	 gen y=ln(Y_hh)
 
	 
	
	*then clcualte efficiency gains
	
	CalculateEfficiency
	 
  restore
end

 
	
	
	
	
	
	
